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It is shown that the formulation of the Einstein equations widely in use in numerical relativ- 
ity, namely, the standard ADM form, as well as some of its variations (including the most recent 
conformally-decomposed version), suffers from a certain but standard type of ill-posedness. Specifi- 
cally, the norm of the solution is not bounded by the norm of the initial data irrespective of the data. 
A long-running numerical experiment is performed as well, showing that the type of ill-posedness 
observed may not be serious in specific practical applications, as is known from many numerical 
p;^ ■ simulations. 
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I. INTRODUCTION 

(N _ 

Taken at face value, the Einstein equations in their original form are second-order equations for the metric compo- 
nents. However, rarely are they used in their original form for numerical applications beyond the harmonic gauge. It 
is much more common to use the 3+1 splitting, introducing the intrinsic and extrinsic curvatures of a spacelike slice as 
£»0 ■ the fundamental variables, in terms of which the Einstein equations become a system of equations that are first-order 
in time and second-order in space [[l]. Some variations are taken, such as decomposing the variables into trace-free 
densities and the traces |0,@|, and, additionally, partially converting the second-spatial derivatives to derivatives of 
first-order variables [Q. 

Alhtough much has been learned in recent years about turning the 3+1 forms further down into a completely 
first-order system |]5|-|l2|, little if anything has been said about these mixed first-in-time, second-in-space forms so 
widely used. 

We investigate here the question of stability of some of these forms under small changes of the initial data, namely, 
the question of well-posedness. We take a particular approach to this question, based on quite standard norms for 
Kr\' periodic solutions. The sense in which the problems that we touch upon are said to be well-posed or stable is the 
following. A system of equations for m variables u(x l ,t), i = 1,2,3, would be stable if the norm of the solution is 
bounded by the norm of the initial data in terms of constants independent of the initial data, namely: 
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IK,*)||<«e bt ||/(.)||, (1) 

where a, b are the same constants for all initial data f(x l ) = u(x l ,0). In the case of linear periodic problems, for 
instance, this means that a, b must be independent of the spectral frequency of the initial data. This definition is 
quite standard 113]. Well-posedness in this or equivalent sense is normally expected of physical systems, as well as of 
successful approximation schemes, including numerical simulations. 

A problem for which the estimate ([!]) does not hold is referred to as ill posed. While showing that the estimate 
holds for all the solutions of a well-posed problem can be accomplished via standard algebraic criteria, showing that 
the estimate does not hold for a certain ill-posed problem may require no more than finding a counterexample. A 
suitable counterexample would consist of a particular sub-family of solutions for which the estimate does not hold. 

By systematically finding "counterexamples" , we are able to observe in Section P that the standard ADM form 0] 
in the two most widely used gauges (geodesic and harmonic slicings) is ill posed in the sense defined above, and that 
its most recent conformally-decomposed version Q is ill posed as well. We may consider this as a contribution to 
the direct analytic investigation of the standard ADM form, which some authors |l(| have termed an outstanding 
problem. 
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This observation might cause a sudden loss of charm to these forms of the Einstein equations. In this respect, it 
is definitely not our intention to question the relevance of ill-posed problems in physics. We refer the reader to page 
230 of H] , where the increasing importance of such problems is rightfully appreciated. After all, the main drawback 
of having to deal with an ill-posed system is surely the current unavailability of relevant results from mathematical 
physics. 



Far from questioning the use of the ADM forms in numerical simulations, in Section III we include a numerical 
simulation of the standard ADM equations with geodesic slicing and periodic boundaries. We show that in spite 
of the ill-posed character, the simulation can be carried out to long times with no signs of instability apparent as 
yet. Furthermore, we refer the reader to Q where numerical simulations with another equally ill-posed version of the 
equations is integrated for long times with no apparent instabilities. We suggest that the ADM forms may retain their 
charm, originating probably from the relatively small number of fundamental variables and apparent compactness 
of the equations, compared to fully first-order forms. Nevertheless, it is likely that numerical simulations of these 
equations might be suject to their inherent instabilites, and an appropriate discretization for long-time evolution may 
not suggest itself in an obvious manner. 

We conclude in Section |v| with some remarks about the reach and relevance of the results. 



II. ILL-POSED ADM FORMS 



Because the reductions of the Einstein equations that we are interested in are rather large systems, we first illustrate 
the procedure that we intend to use in the much simpler case of the 1+1 wave equation. Our procedure is, essentially, 
a rather systematic way of finding counterexamples to ([|) . We take it after an example in page 229 of , credited 
to Hadamard fl5j , (It is hardly contestable that any counterexample will suffice, irrespective of how particularly 
wicked!) 



A. The ill-posed form of the wave equation 

Let's consider the case of the 1+1 wave equation: 

i> = ^xx (2) 

where ' = d/dt and the sublabel x stands for d/dx. Let's reduce the wave equation to first-order-in-time, second- 
order-in-space form by defining a second variable r\ = ip. We obtain thus a system 

f) = ipxx (3a) 
i> = V (3b) 

for a variable u = {-q, tp} with initial values / = {j](x J , 0),tp(x : ' , 0)}. Let's consider, for simplicity, the case of solutions 
of periodicity L. A solution is 

ijj = cos(wi) sin(wx) (4) 
r\ = — uj sin(wi) sm(ujx) (5) 

where u> — 2nn/L and n is an integer. Different initial data are here labeled by different values of u>, and so are 
the solutions arising from them, so we can think that at any given time and position, the value of the solution is a 
function of the initial data through to. Let's calculate the norms: 

1 f L 1 f L 1 



and 



II/(0II = t/ ip(x,0r +r)(x,0ydx = t sm%ujx)dx = - for all 



\u{;t)\\ = \ ip{x,t? +V{x,t) 2 dx 

1 f L 

— —(cos 2 (cut) + uj 2 sin 2 (uit)) / sin 2 (uj x) dx 
L Jo 

= i(cos 2 (^) + to 2 sin 2 (^)) (6) 
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so we have 



||«(-,t)|| = (cosV)+^sin 2 M)||/(.)||. 

Because to 2 has no upper bound in the full spectrum of lu, there are no constants a and b such that 
(cos 2 (wi) + lu 2 sin 2 (a;i)) < ae bt for all initial data (all ui). So we have a subset of initial data for which no esti- 
mate holds, which is sufficient to assert that there is no estimate good for the entire set of data. This form of the 
wave equation is ill posed. 

This is sometimes interpreted as implying that a finite difference scheme for these equations would not progress 
forward in time without instabilities even at short times, on the basis that it is virtually impossible to filter out the 
high frequencies lu [jl3| . 

In the following, we examine two of the most commonly used versions of the 3+1 splitting of the Einstein equations, 
namely the standard ADM form, and the conformally-decomposed version that appeared in [Q . Because the latter is 
more involved, we develop it first in some detail. 

B. Ill-posedness of the conformally-decomposed version of the ADM form 

The system of interest has appeared in jjj, and is based on earlier work in j||. It is a system of 15 equations for 
15 variables (4>, K, 7^ , Ay , T i ), and is referred to as System II in j|, to distinguish it from the standard 3+1 Einstein 
equations jjj, which we display in the next subsection, Eqs. (p7|). These variables are related to the intrinsic metric 
7y and extrinsic curvature Kij as follows: 

(7a) 
(7b) 
(7c) 

(7d) 
(7e) 

where ^ is the inverse of jij . The Einstein evolution equations (^7]) in terms of these variables are equivalent to 
the following: 

d 1 . . 

j^ = - r K (8a) 

j t % = -2aA i3 - (8b) 

j t K = -fW t D 3 a + a (lyA* + ^K 2 ] + l a (p + S) (8c) 

±A\j = {-{D l D J a) TF + a {R% F - Sff)) + a{KA l3 - 2A il A l j ) (8d) 



e 40 


= det( 7y -) 1/3 




Jij 


= e^Ho 




K 


— 7 J^ij 




Aij 






f 







(^' fljl _2^^, m +^ri3 1 ^ . (8e) 



dt 

d 



Here a is the lapse function, (3 l is the shift vector, and f J fc are the connection coefficients of 7^. The superscript 
denotes trace-free part, e.g. Rfj F — i?y —jijR/3. Indices are raised and lowered with 7 y and its inverse. We use the 
shorthand notation 



where £p is the Lie derivative along f3\ Although for the solution to (g) to be a solution to the full set of Einstein 
equations it is necessary that the initial data satisfy the set of four additional constraints, in the following, we do not 
need to consider the four constraint equations explicitly, since they do not affect the well-posedness of the evolution 
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system. We assume that the constraints are imposed on the initial data, thus selecting, from the set of solutions to 
(||), the subset that satisfies the ten Einstein equations. 

The system (j|) is first-order in time and second-order in space, being generically represented in the form 

u = A ij (x k ,t,u)u,ij +B i (x k ,t,u,u ! k)u, i +C{x k ,t 1 u) = P(x k ,t,u,u, k ,d/dx : >)u. (10) 

Any time that an evolution system can be represented in this form, we interpret the right-hand side as an evolution 
operator acting on a solution. In this case, the evolution operator is P, and contains all the terms in the right-hand 
side of (|^) . It is essential to point out that the properties of stability of a system like ( |l0| ) are encoded in the principal 
terms of the operator P, namely, in A lJ (x" 1 , t, u). This means that we can restrict our attention to a system of the 
form 

u = A ij (x i ,t,u)u, i j . (11) 

The lower-order terms that differentiate A 1 ^ (x l ,t,u) from P do not affect the existence of an estimate of the form 
(|l|) (see, for instance, p. 139 of p3| ). If there is such an estimate for (|lT|), then there is one for ([To]) as well, and 
for any other system with a right-hand-side that differs from P only in first and zeroth order terms. If there is no 
estimate for (^TJ), then there is no estimate for ([!(]) either, nor for any system that differs from (^) in first-derivatives 
or undifferentiated terms. 

In the following, we focus on the principal terms of the evolution operator of (^), namely, the exact terms that 
contain the highest derivatives of the dynamical fields in the right-hand-side. In this case, the principal terms are 

(12a) 
(12b) 
(12c) 

-,r'"'"-/., ) ) ■ ( 12d ) 

(12e) 



4> 


= 0, 




= 0, 


K 


= o, 


Aij 


= e-^a ( 


P 





where ' = d/dt. Second-derivatives of the lapse and shift do not contribute to the principal part of the evolution 
operator of (||) if they are assumed to be arbitrarily given as source functions. If they were given dynamically, in 
terms of the metric or extrinsic curvature, then their second-derivatives would contribute to the principal part of the 
evolution operator and should be included (as an example, see the case of the harmonic slicing below). 

For our purposes, it is much simpler to work with system (|l2|) than with (^) without affecting our conclusions. 
The system (|l^) differs from (|h, but only in terms that are of first and zeroth order. Clearly the solutions will be 
different, but not their stability properties. It may not be obvious that the terms that are essential to make ( |12d| ) 
trace- free are of first-order (not second), and therefore do not need to be included in our discussion. The reader may 
as well consider them as included, since their inclusion does not affect our argument in any way. 

We can show that ([l2]) is ill posed by finding one family of solutions for which no estimate of the form ( p]) holds. 
Essentially the same procedure that gave us the result for the wave equation gives the analogous result for ([lj) , and 
consequently, for (^|). However, because the procedure consists of finding explicit sloutions, the results are necessarily 
restricted by the gauge. We will determine ill-posedness in two cases: geodesic and harmonic slicing. For geodesic 
slicing, assume a = 1 and (3 % = 0. Consider the periodic solution 

-log ( 1+ ~cos(k-x) ) , (13a) 



4 V 2 

Jij=Sij, (13b) 
K = 0, (13c) 
(i + cos(k-x)) / 1 



AlJ 4(l + Icos(fc.,)) 3 \ hkj " Z k ' k5 V *' (13d) 
P = 0, (13e) 



where k-x = k l x J Sij, and k l = 2ixn l /L with integer n l . The initial data for this solution are 

(f>(x,0) = -log(l + -cos(k-x)j, (14a) 
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jij(x,0) = 5ij, 
K(x,Q) = 0, 
l y (a:,0) =0, 
f l (x,0) = 0. 



The norm is defined as 



|«(-,t)|| = i / d 3 x(0 2 + ^(7 y ) 2 + ^(r) 2 + x 2 + ^(i y ) ; 



1] I I] 

Let's calculate the norms at the initial time and at time t. We have 



|«(.,0)|| = ||/(.)||= 3- 



1 

1? 



cube 



\ log [ 1 + - cos(fc-x) I I d 6 x, 



and 



IK,t)ll=3 + i£ fee (ilog(l + icos(M)) d 

cos(k-x) + | 
4(1 + ^cos(k-x)) 3 



cube 



X 
2 

d 3 C 



It is our purpose to demonstrate that the norm ||it(-,i)|| can not be bounded independently of 
purpose, we start by factoring out ||/(0||> namely: 



IM,*)|| = ||/(0H i 



2 {k-k) 2 t 2 1 



3 16||/(0||L 3 Lube V(l + icos(fc-a;))3 



cos(fc-a;) + | 



d A x 



Because 1 + 4 cos(fc-a;) < 3/2, we have that 



cos(fc-.x) + | 
(1 + ±cos(/c-:r)) 3 



and because J ube (cos(k-x) + ^) 2 d 3 x = 3L 3 /4 then 



> 



1 

17' 



cube 



cos{k-x) + I 



1 + 5 COs(fc-x) 



cos(fc-a;) + — 



J3 ^ 3 /^ 2 
cTa; > - - 



4 V3 



Using this inequality into ( [i"8| ) we obtain 
Also because 1 + 4 cos(k-x) < 3/2, we have that 



! ^ 2 (fc-fc) 2 i 2 3 (2 



316||/(0|| 4 V3 



> 



log(3/2) V 



which, if plugged into (|Tj), yields 

IK-,*)II> 11/(011 



2(k-k)H 



2/2 



3 /2 



3 16 



log(3/2) ^ 2 4 V3 



Because (k-k) 2 is unbounded, then ||it(-,t)|| increases out of bound at large frequencies u>. This shows that an 
estimate of the type ([!]) does not exist. Therefore the system ( |l2|) is ill posed. The addition of first-derivatives or of 
undifferentiated terms will not turn ( |l2| ) into a well-posed system, from which it follows that (||) is ill posed as well. 
The argument may be impossible to generalize to the case of arbitrary lapse and shift, but presently it suffices to 
make our point . 

Consider now the harmonic slicing a = det 7^ with (3 l = 0. In this case, the principal terms of the evolution 
operator of the system (|J) are not ([[2]), because the second-order derivatives of the lapse are now dynamical and must 
be considered. We actually have a = e 6 ^ and a,ij = 6e 6 ^(0,y +6(j),i4i,j ). Thus the principal terms of the evolution 
of (M) in the harmonic gauge are 



<£ = 0, 
Jij = 0, 

1 



h = e 2 * 



if m 7ii,Im-8U,«-r7«7 ,r ' 



where ' = d/dt. A periodic solution is 



(24a) 
(24b) 
(24c) 

(24d) 
(24e) 



<t>= ol°S 1 + 7;Cos(k-x) 



Jij 



K = 6t- 



cos(k-x)) 



(1 + 5 cos(fc-x)) 



Au = St' 



cos(k-x)) 



+ 5 cos(fc-x)) 



ki k a — k ' k $4 

3 



p = 0. 



(25a) 
(25b) 
(25c) 

(25d) 
(25e) 



With the norm (|15|) and by following very similar calculations as in the case of geodesic slicing, one can readily see 
that the following inequality holds: 



l«(-,*)ll>H/0 



1 + {k-kft 2 



\ 



\ 



59 



log(3/2) V 
2 ) 1 



(26) 



As in the previous cases, we can see that the norm is not bounded because of the presence of u A . This shows that the 
harmonic slicing does not "turn" system (^) into a well-posed form. This appears to contradict standard theorems 
on the well-posedness of the Einstein equations in the harmonic gauge. However, it does not. The difference between 
the use of the harmonic gauge in standard proofs of well-posedness and the treatment in the present paper is that the 
standard proofs are carried out with the Einstein equations either in full first-order form or in the original second-order 
form, whereas here we are considering the mixed case of first-order in time and second-order in space. 



C. Ill-posedness of the standard ADM form 

We now consider the standard ADM equations in the form proposed in Q , which are of wide use in numerical 
applications. The equations are 

7 „- = -2aK l3 + 2D (l f3 }) (27a) 
Kij , = oRi, - 2aK lk Kf + aKK io - D^a + p k D k K l3 - 2K k(i D k /3 j) , (27b) 

where 7^ is the intrinsic metric of the slice at constant t, Kij is the extrinsic curvature of the slice, defined by ( |27t| ), 
a is the lapse function and (3 l is the shift vector. Here we also benefit from restricting attention to the principal terms 
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in the right-hand-side, namely the exact second-derivative terms. Assuming non-dynamical choices of lapse and shift 
(ruling out the harmonic slicing, in particular), the principal terms of the evolution operator of the system are 

jij = 0, (28a) 
k i:j = | (2j kl ji (tJ)k - -r kl j ijt ki - l kL lu, l3 ) ■ (28b) 



Consider the following periodic solution to (|28|) with a — 1 (geodesic slicing) 



With the norm 



we have 



jij = yL + i cos(k-x)J 5i 3 , (29a) 

*« = v+T°X)) {hkj+ s - ,k ' kh <29b) 



||„(.,t)|| s -L j d' x C£(f„f + '£(K„f), (30) 

^ J cube „-„■ 



3-' 



and 

11/(011 = (J J • (32) 
It is very simple to see that the following inequality holds: 

IK^)ll>ll/(-)ll(l + ^(fc-fc) 2 ^) (33) 

Because of the presence of w 4 , we conclude that the norm is not bounded in terms of the initial data and the standard 
ADM equations are ill-posed for a = 1, in the same sense as system (||) is. 

For completeness, we end this section by showing that even the harmonic gauge suffers from the same type of 
ill-posedness. For the harmonic gauge we have a = J det jij > thus DiDja = (a/2)j kl jkl,ij + • • •, so that the principal 
terms of ( p7| ) are 

jij = 0, (34a) 
Kij = | (27 fc S(*j)fc - j kl j ljM " ^J H Jkl,ij) ■ (34b) 

Consider the following periodic solution to (|34"|): 

Jij = ^1 + ^ cos(fc-x)^ Sij, (35a) 
t ( 1 \3 

Kij = — ( 1 + — cos(fc-ir) J cos(k-x) (Akikj + Sijk-k) . (35b) 
With the norm (^0|) it is straightforward to see that the following inequality holds: 

IK,*)II = 11/0)11 (l + ^-fc) 2 * 2 ). (36) 

It is impossible to bound the terms in parenthesis in the right-hand-side by a factor of the form ae bt with a and b 
independent of lu = \Jk-k. Thus an estimate of the form (jl]) does not hold, and the standard ADM equations in the 
harmonic gauge are ill posed. 
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III. NUMERICAL EXPERIMENTS 



We have seen in previous sections that the standard ADM form of the Einstein equations is ill-posed in a certain 
definite sense, and that so are the conformally-decomposed version (|J) and even the 1+1 wave equation in flat space, 
against maybe widespread intuition. Although this may create a difficulty in the stability analysis of these systems, 
we do not think that this necessarily implies an obstruction to numerical integration of any of them, in general. For 
particular applications, it is possible that the bad behavior at high frequencies may be disregarded. On the other 
hand, there may be discretizations suitable for numerical evolution for long enough times, not necessarily arbitrarily 
long. 

In this section we present two numerical experiments with these ill-posed systems. In the first place, we show that 
in the case of the ill-posed form of the 1+1 wave equation there exists a stable discretization. This in principle implies 
that the numerical integration is not hampered in any way. The fact that a stable discretization exists in this case in 
spite of the ill-poscdness is very unusual; it is possible that the 1+1 wave equation be an exception to the rule. 

Secondly, we present a discretization of the standard ADM equations in the geodesic slicing which runs for impres- 
sively long times without any signs of instabilities. 



A. Stable discretization for the ill-posed wave equation 

In this subsection we discretize the ill-posed form of the wave equation, namely (^). This is a standard exercise, 
but a very illuminating one. Our intention is to emphasize that the discretization is stable, which means that the 
discretized equations have no modes that are amplified by the evolution. This is in spite of the ill-posed character of 
the continuous equations that the discretization is modeling. 

We use the leapfrog discretization, which is explicitly given by 

9 At 

'?r 1 = ^ - 2 ^ + + ^ (37a) 

- 2 At??; 1 + (37b) 

where At represents the time step, and h represents the spacing between the grid points, namely h = Ax. The grid 
points are equally spaced in the interval —1 < x < 1, Xi = — 1 + (i — l)h for i = 1 . . . N and h — 2/(N — 1). 
In order to carry out a stability analysis [ fl6| , we use the ansatz 

V'J = Ce llkh Vo (38) 
W = Ce ijkh ipo (39) 



where i = y— T. Substituting this ansatz into ( 37a ), we obtain 



(e - 1) ~4^(cos(fc/i) r,„ 

-2 At£ (£ 2 -l) )\^o)~ [ ] 

For a non-trivial solution we need the determinant of the system to be zero, namely 

{( 2 - if + 2a 2 {l -cos{kh))( 2 = (41) 



where a = The solution is 



e 2 = -(c-l)±V(c-l) 2 -l (42) 



with c = a 2 (l — cos(kh)). The discretization is stable if |£| < 1. A sufficient condition for this to happen in our case 
is that c < 2, where we have = 1, i.e. the discretization is not only stable but unimodular. The requirement for 
c < 2 is that a < 1. This means that for our discretization to be stable and unimodular we only need to take the 
time step smaller than one- half the grid size h. 



We implemented a straightforward Fortran code based upon the discretization (37a), with periodic boundary 
conditions in the domain — 1 < x < 1 . As expected, the code reproduces very well the evolution of a pulse of compact 
support of the form 

— ((t-vx f - w 2 ) for \t-vx\<w, ^g-j 
otherwise, 
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where the evaluation of the condition is to be understood modulo 2 (the periodicity of the grid) . The pulse shown in 
Fig. |l| has amplitude A = 1 and width w = 1, and it is propagating to the right with speed v = 1. We can see no sign 
of damping or amplitude growth, and very little distortion, even after 100 crossing times. Shown in the figure are the 
initial time, t = 0.0 (solid line) and the final pulse, t = 200.0, at a resolution of N — 200 points (dotted line) and 
N = 400 points (dashed line) respectively. Note how the distortion of the pulse decreases with increasing resolution 
as expected from the second-order discretization (37a). 



B. Numerical stability of the ADM equations 

The ADM equations in the form introduced by York Q, and variations of the same, have been used extensively in 
numerical work [^|-| 17-2^], to cite a few. For completeness, however, we present here an alternate straightforward 
numerical implementation of the ADM equations, which illustrates that one can indeed integrate the nonlinear system 
(not just its linearized counterpart, and not just its principal part). Starting from J27j), we restrict our attention to 



the case of geodesic slicing, where a = 1, /3 ! = 0. The calculation of the Ricci tensor in (27b) is the most involved 
and error-prone task, but from the standard definition ]2l| l 

p. . pk I pfe pro pni I pn jnm (A-Ast ) 

^ij 1 ij,k ' ik.j 1 im L jn ' 1 ni mj ' \ <*) 

= \l M + + lli,j) , (44b) 

it follows that its computation can be organized by collecting terms which contain second derivatives, terms with first 
derivatives of the metric and its inverse, and terms with only connection components, 

R V = -jl kl ( 2 7l(i,j)k - Hj,kl ~ lkl,ij) 

+ ,fc {-Hj,i + 7iM + lu,j) - -jl ,i {-Jkj,i + 7ji,k + llk,j) 

-pn pm _r_ pn pm / az\ 

im jn ' ni mj' V / 

We first compute the inverse metric and the second spatial derivatives of the metric, from which we obtain the first 
term, then calculate the first derivatives of the metric (and of the inve rse metric) , which yields the second and third 



term. The final step is to calculate the connection coefficients (44b) and the last two terms in (|45|). The terms 
involving Kij in the right hand side of ((2^) are expressed in terms of the (downstairs) extrinsic curvature and the 
inverse metric, 

KKij - 2K ik K) = 7 ' m (K lm Kij - 2K H K mj ) . (46) 
We discretize (^7|) using ( 44b )-(f45|), taking the variables 7jj and Kij as given on a rectangular, equally spaced grid 



of size N 3 , with resolution h = 2ir/N and covering the domain < x < 2tt, < y < 2it and < z < 2ir. We 
label the grid points by = k h, yi = I h, z m = mh, with k,l,m = 1 ... A, and the time levels as t — to + nAt, for 
n = l...A t . To distinguish tensor indices from grid indices, we use the notation Jij[ki m ] = lij{ x k,yi, z m ,t n ). We 
use symmetry where appropriate, storing only the relevant components of 7y, Kij, 7 y , 7y.fc, 7 y ,fc, "fij,ki and Rij. For 
instance, in the calculation of the second and third terms of ([45|), the r y' l \k are stored on the space allocated later to 
therf,.. 

We compute first spatial derivatives jij.k and 7 y ,fc, and second spatial derivatives jij.ki with centered, second-order 
accurate finite differences on grid points, e.g. 

= ^ (V'G+ij.fc] " K-i^k}) (47a) 
t«*Kj,k\ = 7^2 - 2 Wj,k] + W-ij,k]) (47b) 

^y[lj,k\ = ^2 &\bi d+ i, k ] - W-ij+i,k] ~ i>li + i,j-i,k\ + ^-i,j-i,k}) (47c) 

We perform the grid indices operations identifying the points xq = xn and xn+i = xi, which enforces periodic 
boundary conditions in all coordinates. 

The time integration scheme we use is the so-called iterative Crank-Nicholson (ICN) method |l7|JT^ fl, with three 
iterations (i.e. one predictor step "forward in time", plus two correction steps). Considering ( p7| ) as equations of the 
form u ~ Pu, where P is an operator acting on« = {lij, Kij), the time integration algorithm is given by 







< +1 =u n + AtPu n (48a) 

u 2 +k = \{< +1 +u n ) (48b) 
u^ +l =u n + AtPuJJ** (48c) 

(48d) 



2 

u n+1 =u™ + AtPu3 + ^ (48e) 



where quantities with sub-indices are intermediate values which do not require additional storage. 
To test the algorithm, we give as initial data the following, evaluated at t = 

jij = [A (miirij — riiTij) + B (rriinj + riirrij)] sin(fc-x — ujt) + Sij (49a) 
Kij = [A (m imj - n inj ) + B (m inj + n imj )] co S (k- x - ort)| (49b) 

where m l and n l are unit vectors orthogonal to the propagation vector k l , i.e. m-m = 1, n-n = 1, m-n = 0, m-k = 
and n-k = 0. For k — (ml,m2,m3) with ml,m2,m3 integers, Eq. ( |49] ) is a periodic solution of the linear system 
obtained by setting 7^' = S 1 ^ in (|27j). During the integration, we monitor the Hamiltonian constraint, given by 

H = R + K 2 - K tJ K l] (50) 



and, as a measure of the stability of the algorithm, the norm given in Eq. (30). We have followed the evolution of a 
pulse with A = B = 10~ 6 , k = (1,2,1), lo = ^6, taking m = (-7,4, l)/v4i6 and n = (-1,-1,3)^6/11. The solution 
is periodic in time as well, with period T — 27r/v / 6, which determines the natural time scale for the test problem. The 
equations can be integrated on a grid of (48) 3 points for hundreds of crossing times, without any signs of instability, 
as evident from Fig. I, which shows ||u(-,*)||/||/(-)||, the norm of the numerical solution as a function of time, relative 
to the initial norm. The relative change in the norm, with respect to the initial norm, is shown as well, and is below 
one percent at 100 crossing times and « 19,000 iterations. 

We have plotted in Fig. || the maximum absolute value over the grid of the Hamiltonian constraint, Eq. (|5(]) as a 
function of time, and this quantity remains bounded throughout the evolution, for upwards of 19,000 time steps. 

We have not carried out the Von Neumann stability analysis for the principal part of the ADM equations, as we 



did for the simpler case of the wave equation in subsec. Ill A. since it would be rather involved and well beyond the 



scope of this work. To our knowledge, this type of analysis has not been made for the other systems mentioned here, 
either. 



IV. REMARKS 



In the first place, it is essential to emphasize that we have kept ourselves within the context of problems with periodic 
boundaries, for considerations of analytical and numerical stability alike. In the case of the Einstein equations, more 
often than not instabilities develop in the course of a numerical simulation, but the origin of the instabilites is not 
known, being attributed intermittently to either the equations themselves or the particular boundary conditions 
being used. In the line of a number of authors (most recently, ]l7[-p0|) who have used the standard ADM equations 
in exact (nonlinear) form and their conformally-decomposed versions in numerical relativity, including evolution in 
the strong-field regime for blackhole mergers, we have presented a long-running stable code for the exact standard 
ADM equations up to (un)stable boundaries. With such a code, which runs for sufficiently long times with periodic 
boundaries even if the equations themselves are ill posed, the origi n of instabilitites, if they occur, can be shifted 
to the boundaries. A discretization such as that used in subsection IIIB might be used to identify stable boundary 
conditions other than periodic, or it might be used for matching to a stable exterior characteristic code j22]-|24]|, in 
which case the boundary values are provided by the exterior code and are, in principle, consistent with the interior 
solution. 

We have not been involved with constraint propagation in this work, but in this respect it is worth emphasizing that 
not only does the code run for long times for the standard ADM equations, but it does so preserving the hamiltonian 
constraint as well, in spite of the ill-posedness of the evolution equations. 
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Secondly, we have shown here that, from the analytical point of view, there is no advantage to the conformally- 
decomposed equations (J] with respect to the standard ADM form 0]. Although a striking difference between the 
numerical integration of the two systems has been reported Qj, the origin of the difference in numerical behavior 
must lie in some factor other than well-posedness (sec |llj for a discussion of other factors). It might be thought 
that both systems may have different properties when turned to full first-order form, and that this might explain the 
difference in numerical behavior. So far this is an open question. In this respect, it has been shown that reducing 
the conformally-decomposed system (||) down to full first-order form by defining the spatial derivatives of the metric 
as new first-order variables does not automatically make it well-posed, unless the lapse function is densitized and the 
constraints are combined in specific ways with the evolution equations p"l| . Nevertheless, in reducing to first-order 
form, there is plenty of freedom in the choice of first-order variables. It is not clear at this time whether or not there 
exists a choice of first-order variables for the reduction of (|8|) which will turn the system into a well-posed one without 
densitizing the lapse or combining with the constraints. The same might be said about the standard ADM form, 
Eqs. @. 
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FIG. 1. Evolution of a pulse traveling to the right at speed v = 1. Shown in this figure are the initial (t = 0.0) and the final 
pulse (t = 200.0). The later time is seen at the two grid resolutions, N = 200 (dotted line) and N = 400 (dashed line). 
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FIG. 2. Evolution of a linearized, periodic ADM pulse traveling in the direction k — (1,2, 1). Shown in this figure are the 
scaled norm N/N = ||«(-,*)||/||/(-)ll ( solid line ). and its relative change N/N -1 = ||«(-,*)||/||/(-)ll ~ 1 - (dotted line). Note 
that the relative change in the norm is below one percent at t = 200-7T, that is 100 crossing times and w 19, 000 iterations. 
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FIG. 3. The maximum absolute value over the grid of the Hamiltonian constraint as a function of time for initial data 
corresponding to a linearized, periodic pulse. 
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